Stochastic statistical theory of nucleation and evolution of nano-sized precipitates in 
alloys with application to precipitation of copper in iron 
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The consistent and computationally efficient stochastic statistical approach (SSA) is suggested 
to study kinetics of nucleation and evolution of nano-sized precipitates in alloys. An important 
parameter of the theory is the size of locally equilibrated regions at the nucleation stage which is es- 
timated using the "maximum thermodynamic gain" principle suggested. For several realistic models 
of iron-copper alloys studied, the results of the SSA-based simulations of precipitation kinetics agree 
well with the kinetic Monte Carlo simulation results for all main characteristics of microstructure. 
The approach developed is also used to study kinetics of nucleation and changes in microstructural 
evolution under variations of temperature or concentration. 
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INTRODUCTION 



Studies of microstructural evolution in phase- 
separating alloys attract interest from both fundamen- 
tal and applied points of view. ^From the fundamen- 
tal side, elucidation of microscopic mechanisms for the 
formation and evolution of embryos of new phases aris- 
ing under the first-order phase transitions is one of the 
principal problems in statistical physics being not well 
understood as yet 1-4]. From the applied side, under- 
standing factors that determine different characteristics 
of microstructure formed under precipitation is impor- 
tant to control these characteristics, particularly for al- 
loys with nano-sized precipitates which attract recently 
much attention in connection with industrial applications 



Presently, theoretical studies of the precipitation kinet- 
ics employ usually either the phase-field method (PFM) 
[j-flij or Monte Carlo modeling (j. |]J5r{l8| . However, em- 
ploying the phase-field methods to describe nucleation 
and evolution of nano-sized precipitates can be mislead- 
ing for at least three reasons. First, the "continuous" 
approximation used in the PFM disregards the discrete 
lattice effects which should be important at first stages 
of nucleation when typical precipitate sizes are few lat- 
tice constants. Second, at the concentration and tem- 
perature values typical for applications, the mean-field- 
type CALPHAD expressions for thermodynamic poten- 
tials usually employed in the PFM studies [§-[l3| will be 
shown to strongly distort the position of spinodals, thus 
using these expressions can drastically distort the type 
of microstructural evolution. Third, treatment of fluc- 
tuative terms (being crucial for the nucleation stage) in 
the "stochastic PFM" versions used in applications until 
now [U[I(J seems to be arbitrary and inconsistent (l9ll20j . 
while the adequate description of these terms determines 
all main characteristics of microstructure. 

Therefore, the only reliable source of theoretical in- 
formation about nucleation and evolution of nano-sized 
precipitates now is Monte Carlo modeling, in particular, 



the kinetic Monte Carlo approach (KMCA) developed 
in Refs. [J, [l5rlT8j . However, present versions of the 
KMCA are time-consuming, which may partly explain 
a relatively few number of applications of this method 
to concrete systems [l7l Il8j|. Moreover, the lattice mis- 
fit and elastic strain effects important for many phase- 
separating alloys can not be easily taken into account in 
the KMCA, while it makes no problem for statistical ap- 
proaches [21(. Finally, in the KMCA it is often difficult 
to follow the dependence of various characteristics of evo- 
lution on thermodynamic and microscopic parameters of 
an alloy, such as the composition, temperature, differ- 
ent interaction constants, etc, while it is usually much 
simpler for statistical methods based on some analytical 
equations. Therefore, the development of a consistent 
statistical theory which takes into account all achieve- 
ments of the KMCA seems to be very important for a 
more deep understanding of the phase separation kinet- 
ics. 

Recently Sroev et al. [2(| (SPV) presented an attempt 
to develop such a theory using the stochastic statisti- 
cal approach (SSA) described below in Sec. Ill B. To 
illustrate the main ideas of this approach, SPV used 
only simplest methods and models, such as the mean- 
field approximation (MFA), continuous approximations, 
the direct-atomic-exchange rather than the vacancy ex- 
change kinetic model, oversimplified interaction models, 
etc, while no attempts of quantitative treatments for re- 
alistic alloy models have been made. 

The main aim of the present work is to raise the accu- 
racy and the predictive power of the SSA in describing 
the precipitation kinetics for realistic alloy models up to 
the level comparable to that of the KMCA. To this end, 
we perform detailed studies of nucleation and evolution 
of nano-sized precipitates in Fe-Cu alloys using both the 
KMCA [H [H and the SSA. This requires many refine- 
ments of simple models used by SPV. We have to consider 
the realistic vacancy-mediated exchange kinetics rather 
than the simplified direct-atomic-exchange model; to use 
the quantitative, cluster statistical methods rather than 
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the simple MFA; to allow for strong concentration and 
temperature dependences of generalized mobilities in the 
resulting kinetic equations, etc. All these refinements 
are made in the present work. We also introduce the 
"maximum thermodynamic gain" principle to determine 
the key kinetic parameter of the SSA, the characteristic 
length of local equilibrium in the course of the nucleation 
process. 

In Sec. II we discuss the methodological problems men- 
tioned above: generalizations of our statistical approach 
to the vacancy-mediated kinetics case; employing cluster 
methods for both thermodynamic and kinetic statistical 
calculations; elaboration of effective methods for calcu- 
lations of effective mobilities which enter resulting ki- 
netic equations, etc. Here we also generalize the earlier- 
suggested "equivalence theorem" |22j | which enables us 
to reduce the vacancy-mediated kinetics to that for some 
equivalent direct-exchange models. In Sec. Ill we re- 
mind basic ideas of the classical theory of nucleation and 
present the main equations of the SSA. In Sec. IV we dis- 
cuss the models and the methods of simulations used and 
we describe the "maximum thermodynamic gain" prin- 
ciple suggested to estimate the local equilibrium length 
mentioned above. In Sec. V we discuss the features of mi- 
crostructural evolution observed in the KMCA and SSA 
simulations for various alloy states considered. Here we 
also use the SSA to study kinetics of nucleation and in- 
fluence of variations of temperature or concentration on 
evolution of microstructure. Our main conclusions are 
summarized in Sec. VI. 



II. QUASI-EQUILIBRIUM KINETIC 
EQUATION FOR VACANCY-MEDIATED 
KINETICS IN SUBSTITUTIONAL ALLOY 

A. General equations for mean occupations of 
lattice sites 

In this section we derive kinetic equations for mean 
occupations of lattice sites disregarding fluctuations of 
atomic fluxes which lead to local violations of the sec- 
ond law of thermodynamics. These "quasi-equilibrium 
kinetic equations" (QKE) differ from the stochastic ki- 
netic equations discussed below (which take into account 
such fluctuations) and generalize those used by SPV [20[ 
for simplified direct-atomic-exchange models to more re- 
alistic, vacancy-mediated-exchange (VME) models. Here 
we also generalize a similar treatment of the VME kinet- 
ics made by Bclashchenko and Vaks 221 (BV) to more 
realistic VME models used in Refs. [l7| and below. 

We consider a substitutional alloy with (to + 1) com- 
ponents p' which includes atoms of m different species 
p = Pi, P2, . . . p m and vacancies v: p' = {p, v} . The dis- 
tributions of atoms on the lattice sites i are described by 
the different occupation number sets {nf } where the op- 
erator n% is 1 when the site i is occupied by a p'-species 
component and otherwise. For each i these operators 



obey the identity J2 P ' n i = 1; so OIU y m 01 them are 
independent. It is convenient to mark the independent 
operators with greek letters p or a : (nf )mdc P = n?, 

while the rest operator denoted as is expressed via 
p 

n i- 

n? = l-5>?- (1) 

p 

For dilute alloys, it is convenient to put "ft" in (TT]) to be 
the host component, e.g., h—Fe for FeCuv alloys. 

In terms of all operators the configurational Hamil- 
tonian H' (for simplicity supposed to be pairwise) can be 
written as 

H ' = \ E y^<^- (2) 

p'q' ,ij 

After elimination of operators according to (p}, the 
Hamiltonian ([2]) takes the form 

H = E + £ <Pp n i + H int 
pi 

Hint = £ (3) 
pcr,i>j 

which includes only independent n^, while constants E , 
ip p , and "configurational interactions" are linearly 

expressed via the interactions V? q in (|2[). in particular: 

v u = (y pa - y ph - y ha + v hh h ■ ( 4 ) 

The fundamental master equation for the probability P 
of finding the occupation number set {rof } = £ is [2X| ] : 

dP{0/dt = J2[W(^,V)P(V) - W( V ,OP(0] ee SP (5) 
v 

where W(£-, rj) is the rj — ¥ £ transition probability per 
unit time. Adopting for probabilities W the conven- 
tional "transition state" model [3, Il5l - fl7j , we can express 
the transfer matrix S in ([SJ in terms of the probability 
of an elementary inter-site exchange ("jump") pi ^ vj 
between neighboring sites i and j: 

W?; = nfnj^ pv exp[-f3(E^ p VJ - E% <yj )] (6) 

where uj pv is the attempt frequency, j3 = 1/T is the 
reciprocal temperature, Ep f V j is the saddle point energy, 

and E™ v j is the initial (before the jump) configurational 
energy of a jumping atom p and a vacancy. The saddle 
point energy E^f v j, generally, depends on the atomic 
configuration near the bond ii (which is neglected in 
simplified kinetic models [3, [H, [22| ) . We will describe 
this dependence by the model used in Refs. [l6[ and 
[ItJ supposing the saddle-point energy to depend only 
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on occupations of lattice sites I nearest to the center of 
bond ij (these sites will be denoted l% n ): 



E, 



SP 



(7) 



Here is the saddle point energy for the pure host 
metal, while the operator A p - describes changes in this 
energy due to a possible presence of minority atoms near 
the bond: 



A p - 



E 

P , z=4 J , 



(8) 



where 



is the total number of nearest lattice sites I 



for each bond (being z nn = 6 for the BCC lattice), and 



As discussed in detail in [21], for the usual conditions 
of phase transitions corresponding to absence of external 
fluxes of particles or energy (that is, when the alloy is a 
"closed" but not an "open" statistical system), the dis- 
tribution function P(£) = P{n^} in ([5]) can be written 
as 



P{<} = exp 



(9) 



Here parameters A^ (being, generally, both time- and 
space-dependent) can be called "site chemical potentials" 
for p-species atoms; Hint is the same as in ([3]); and 
the generalized grand canonical potential f2 is determined 
by normalization. The analogous equation (10) of BV 
[22| (who treated both closed and open systems) differs 
from ([9]) with replacing the interaction hamiltonian by 
a more general "quasi- interaction" operator Q. For the 
closed systems discussed in this work, Q = H int , and 
this relation greatly simplifies kinetic equations discussed 
below with respect to those of BV. 

Multiplying Eq. ([5]) by operators n p and summing over 
all configurational states, i.e. over all number sets {n? }, 
we obtain the set of equations for mean occupations of 
sites ("local concentrations") cf = (n?): 



dc^/dt = MS) 



(10) 



where ((...)) — Tr{(...)P} means averaging over the dis- 
tribution P, for example: 



4 = «> = E C * P K}- 



(ii) 



After a number of manipulations described by BV, 
Eqs. (fT0| can be transformed into the QKE for mean 
occupations cf . These equations are similar to Eqs. (19) 
and (42) of BV but include generalizations and simplifi- 



cations mentioned above: 

def/dt = X><" 6 §(W-W) (12) 

dej/dt = [^^v^+E^vfef^f) 

-{*-»• J'}]- (13) 

Here and below, greek indices a, /3 . . . correspond to the 
minority atoms; symbol "v" is used for vacancies; in- 
dex nn means "nearest neighbors" , and symbol u jn n (i)" 
means summation over sites j being nearest neighbors 
of site i. Term 7 pv in (fT3|) (where p is a or ft, i. e. 
a minority or host atom) is the effective exchange rate 
p f± v for a pure host metal. This term can be written 
in the form similar to Eq. (j6|): 



7 P v 



exp(-/?££Q 



(14) 



wher 



is the same as in ([6]), while is the effec- 
tive activation energy which is expressed via the saddle 

© as 



point energies E£ in 
follows: 



and interactions V- 



p q 



E pv 

etc 
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(V ph + v vh ) + v h 



(15) 



Note that for minority atoms p =a, expressions (|15p dif- 
fer from analogous activation energies P P J MC used in 
the KMCA and given by Eq. (2.5) in Ref. 



pQV . TPQTV 



Kl = Kluc- (16) 



The difference arises because in the statistically averaged 
Eqs. (fT2|) and (fT3|) . the probability ([6]) is averaged over 
distribution ([9]) , and for the inter-site exchange ai ^ vj , 
it leads to the additional Gibbs factor exp(— /3v£j v ) in 
the averaged probability. 

Quantities in (| 13j) (for brevity to be called "cor- 
relators") are certain averages of site occupations which 
describe influence of minority atoms in vicinity of the 
bond ij on the pi^vj jump probability: 



al 



,1=11; 



(17) 

where is the same as in (|SJ), while quantities uf t (to 
be called "kinetic interactions" as they affect only effec- 
tive jump probabilities but not thermodynamic proper- 
ties) are related to the interactions V? q in ([2]) as fol- 
lows: 

«s = n h - v$ h . (is) 

Finally, quantities £J and rjf in (fl3|) and (fT2ll can be 
called "site thermodynamic activities" for vacancies and 
a -species atoms, respectively, as they are related to site 

as: 



chemical potentials Af in 
^ = exp (^A?); 



r,f = exp(/3Af), 



(19) 



that is, analogously to the relations between conventional 
thermodynamic activities and chemical potentials. 
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B. Calculations of site chemical potentials A? and 



correlators bf- 



To solve QKE (|13[) . we need explicit expressions for 
site chemical potentials = A^(cj) determined by Eqs. 
(|TT|) . and for correlators = fey(cfc) determined by 
Eqs. ([T7| . To find these expressions, we should use some 
approximate method of statistical physics, such as the 
MFA or cluster methods [2l[. As discussed in detail in 
[241 125I ] and below, at realistic values of interactions v"? 
that significantly exceed temperature T (which is typi- 
cal, in particular, for the iron-copper alloys under consid- 
eration) , employing the MFA leads to great errors in cal- 
culations of thermodynamic potentials which exclude any 
realistic description. At the same time, the pair cluster 
approximation (PCA) usually combines simplicity of cal- 
culations with a rather high accuracy, particularly at low 
c and T under consideration, see, e. g. [241426^ . Thus, 
for site the chemical potentials Af we use their PCA- 
expressions which for a binary alloy ABv with minority 
atoms A, host atoms B, and a realistically small concen- 
tration of vacancies: cj <C 1, are given by Eqs. (39) of 
BV: 



X? =T 



]n(ci/<£) + X>(l 



Here Aj = A", < 



ln(cr/cf) 



E 



(20) 

_ (21) 

i\ <k = 4, 4 = (l-ci-cj) ~ (l-a), 

while the function gtj or gjj is expressed via the Mayer 
function /y = [exp (-fay)-!] or = [exp {-^vJ A )~ 



1] for the potential Vi 



— „,AA 



or v 



vA 



defined in 



vA 



vA 



2^ 



AB 



V; 



BB 



AB 



BB 



(22) 



as follows: 



9ij = 2fij/[Rij + 1 + fij{(H + Cj)] 
9 r ii = 1fl j f[Ri j + l + fti(c i -c i )] 

Rij = {[1 + + c*)/y] 2 - AciCjfijifij + 1)} 



1/2 



(23) 



Let us also present the PCA expression for the free 
energy F of a binary alloy [2l[ which is used below in Sec. 
IV B for discussions of thermodynamics of precipitation: 

F = ]T T In c\ ~ t\ ln(l - !/,,<;<■, ; + £ A 4 c 4 (24) 



If temperature T, 
T » wy, the PCA ex- 



where gy is the same as in 
much exceeds interactions Uy 

pressions (|20|) and (|24)l transform into the MFA ones [21 1 . 
However, as mentioned, in situations of practical interest 
we usually have an opposite inequality: T <C My, and 
employing the MFA is misleading. 



Let us now discuss calculations of correlators bf 4 
in p7|) . For simplicity, we first consider the case 
of configuration-independent saddle-point energies when 
differences A£ in Eqs. ((8]) and dTTl) vanish and correla- 
tors &?■ = 6y do not depend on the kind of a jumping 
atom p. Using identities 

{nff=nf, nfn^ = nfS aP , nf nf = 0, 

exp (sn?) = l +</(*), f(x) = (e x -l), (25) 

we can rewrite Eq. (fT7|) for this case as follows: 



;=i a 



fet 



E E E (W-^/C-C ( 2 6) 

fc=0 h^...lk oti...a>k 



where quantity /" is defined as 



(27) 



with /(x) from (f2"5")l . while k t in (|26p is the total num- 
ber of sites with nonzero values of potentials (ug + 
For example, for the nearest-neighbor or second-neighbor 
interaction models in the BCC lattice, we have: k t = 14, 
or: kt = 20. 

To find averages in (|2l?)) over distribution (JSJ) , we should 
again employ some approximate method of calculations, 
such as the MFA, PCA or higher-order cluster approx- 
imations [2l], (26|. However, for the most of systems of 
practical interest, in particular, for the iron-copper al- 
loys discussed below, we can apparently use in (1261) the 
simple MFA replacing each operator nf by its average 



value c? 



1 1 seems to be justified as the functions 



ff in Eqs. (f2"6"|) and (|2T|) for such systems are typically 
rather large (for example, f{j3u\) > 5 for the systems de- 
scribed by Table III below). Thus the main contributions 
to sum (126|) come from averages of products of many dif- 
ferent operators nf which correspond to well-separated 
and weakly correlated sites I. In particular, for the BCC 
lattice, these products (even for the nearest-neighbor in- 
teraction model) include terms with the neighbors from 
first to tenth, most often third and fourth. Correlations 
of occupations for the so distant lattice sites should be 
typically small, thus using the MFA that neglects such 
correlations should be adequate. 

In the MFA, Eq. fl26J) takes the form 



c h c h 



E E E 

k=0 l\=jt...l m =jti,j ai...O| 



f<*k 

■Jh. ■ 



(28) 

To further simplify this expression, we can take into ac- 
count that the space variations of local concentration cf 
arising in the course of alloy decomposition are typically 
rather smooth, particularly at the nucleation stage (for 
which an adequate description of kinetic coefficients that 
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include correlators 6y is most important), see, e. g. Figs. 
|2"0"1I2"21 below. Thus, the c? values for sites I adjacent to 
bond ij that enter Eq. ([28f are usually close to both the 



(29) 



cf and cf, as well as to their average 



c«)/2. 



Therefore, to avoid unnecessary complications of compu- 
tations, we can approximate each cf in Eq. (|28|) by the 
average (|29l) . It enables us to write the correlator 
in the simple analytic form: 



w n a 



(30) 



For example, for a binary alloy ABv in the BCC lattice 
described by the second-neighbor interaction model with 
two kinetic interaction constants, u\ and 112, we obtain: 



h h 
c, c 



6 



X 



' Ctf/GSuOHl + cy/^Ua)] 6 



(31) 



where index a =A in and is omitted for brevity. 

When differences A? in Eqs. ([8|) and (fT7j) are nonzero, 
the correlator &?■ in Eq. (TiTl) can be calculated by the 
same way as bij in Eqs. ([26l) - ([3"Tj) . The difference arises 
only for sites I = l% n adjacent to bond ij for which fac- 
tors /" defined by Eq. ([27jl are now replaced by analo- 
gous factors defined as: 



rap 

Ja 



(32) 



In particular, for the BCC binary alloy ABv with the 
second- neighbor interaction, we obtain instead of (|31|) : 



b^= ^[1 + 5,,/aTx 



where /^ p = 



f[P(ui+U2 



(33) 



A?)], and Ap = ( £ p - e p A ) 



C. Equivalence of precipitation kinetics for the 
vacancy-mediated exchange models to that for 
certain direct exchange models 

In this section we show that the VME kinetics de- 
scribed by Eqs. (fT2")l and (TT3"|) can usually be described 
in terms of certain equivalent direct-atomic-exchangc 
(DAE) models. It will generalize the analogous "equiva- 
lence theorem" derived by BV. 

First we note that the vacancy activity Q = exp(/3Ay) 
in Eqs. (fT2"j) and ([T3"|) is proportional to the vacancy 
concentration cj. It is illustrated by Eqs. (|2"Tj) and is 
actually a general relation of thermodynamics of dilute 
solutions. Thus time derivatives of mean occupations are 
proportional to the local vacancy concentration cj or cj, 
which is natural for the vacancy-mediated kinetics. As 



cj in real alloys is very small, this implies that the re- 
laxation times of atomic distribution {cf} are by a fac- 
tor 1 /cj greater than the time of relaxation of vacancies 
at the given {cf } to their "quasi-equilibrium" distribu- 
tion cj{cf} for which the right-hand side of Eq. (p~3|) 
vanishes. Therefore, discarding small corrections of the 
relative order cj -C 1, we can rewrite Eq. ([T3"l) as follows: 







JnnW 



hwb% 



E^g^-O^i}]- ( 34 ) 



which can be called "the adiabaticity equation" for the 
vacancy activity Solving this equation we can, in 
principle, express Q via c" . Then substitution of these 
£^(c") into Eq. (fT2|) yields the QKE for some equivalent 
DAE model. 

To illustrate this approach, we first consider the VME 
models with configuration-independent saddle-point en- 
ergies. For such models, parameters Ap in (|8|) are zero, 
correlators 6 P = bij do not depend on the kind p of 
jumping atom, and the adiabaticity equation (|34|) takes 
the simple form: 



E b ^ 



7hv+52l*vri?)/&-{i->j}\ =0 (35) 



If we denote the first term in square brackets (|35j) as 
1/vii then the difference in these brackets takes the form 
vf l — Thus the solution of Eqs. ([331) is provided 

by Vi being a constant independent of the site number i 
(though possibly depending on time, as well as on tem- 
perature and other external parameters): 



lavV; 



') = v(t). 



(36) 



Relation ([56]) determines the above-mentioned "quasi- 
equilibrium" vacancy distribution c^{cf} which adiabat- 
ically fast follows the atomic distribution {cf}. Substi- 
tuting it into Eq. (|12j) we obtain the explicit kinetic 
equation for atomic distributions {cf } for which influ- 
ence of vacancies is characterized by a single parameter 
v(t) being a "spatially self-averaged" quantity: 



dcf/dt = E 



Jnn(») 



E 7*v7^v (rjj 




a P 

Vi - V i Vj 



(37) 



The last term of this equation (missed in the analogous 
Eq. (46) of BV) is present only for many-component 
alloys with two or more species of minority atoms. Eqs. 
(37) can also be rewritten in the form used for DAE mod- 

dcf/dt = E M% h 2sinh[P(\? - \f )/2] 

+ E M^2sinh[^+Af-Af-Af)/2] (38) 
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where generalized mobilities M? q which describe inter- 
site exchanges a^p and are given by the following 
expressions: 



In the case (b), we can rewrite Eq. (|42l) as: 



MS h = 7a v7fcvK*) K exp + Xf)/2] 



(39) 



= 7 av7/3v^(i) h 3 exp [f3(\f + A" + Af + Ap/2]. (40) 



Comparing these expressions to Eq. (32) of BV for mo- 
bilities M^ q in an alloy with the nearest-neighbor direct- 
exchange rates 7^ q =7pq, we see that Eqs. ([55]) and (T4T)]) 
correspond to a DAE model with the following effective 
direct exchange rates: 



llh = 7av7hv V(t)\ 7^ = 7av7/3v v(t) 



(41) 



Note that the effective DAE rates ([41]) are by a factor cv 
smaller than the vacancy exchange rates 7 pv . 

Let us now consider more realistic VME models with 
the configuration-dependent saddle-point energies when 
correlators in ([17]) for different p are different. For 
such models, the basic adiabaticity equation (|3~4"|) for va- 
cancy activities generally, can not be solved analyt- 
ically, thus either numerical or some approximate ana- 
lytical methods should be used. Let us discuss two such 
approximate methods employed below. For a binary al- 
loy, we can rewrite Eq. (|34[) in the form: 



E h % s^Thv [(i + '/>'•,.) /g • '/,'•,) /g 



= 
(42) 

where ?7i=7?f =exp (/3A»), and r„ is 7av&g-/7fev&y- E Q- 
(1421) can be approximately solved if products r^^f obey 
either of two inequalities: 

(a) runj = exp(/3Ai)7 QV 6 4 Q j /7/ lv &^ < 1; 

(b) mrij = exp(/3A J )7av&y /7ftv&' i i » 1. (43) 

In the case (a), the second terms in round brackets in 
(f4"2"j) are just small corrections to the first ones. In the 
zcroth approximation they can be neglected, thus the 
zero-order solution of Eq. (|4"2l is: 



(i) = u(t) lhv 



(44) 



where the constant factor 7^ v is introduced so that the 
function v(t) is analogous to that in (|36[) . Substituting 
(0 into (fI2|). we again obtain Eqs. ([37]) or (J35J for a 
binary alloy: 

dc z /dt= E My2sinh[/3(A 3 - AO/2]. (45) 

Here index "ah" at the effective mobility Mf- 1 = M id 
is omitted for brevity, the expression for this mobility is 
similar to that in Eq. (|39l) : 



M«(a) = 7a" 6« ex P ^( A * + A i)/ 2 ]' 



(46) 



Jn„(i) 



LQ £V £V ~, 



= 



(47) 



where the second terms in round brackets are again small 
corrections to the first ones. Thus the zero-order solution 



of this equation can be written as: ^ 



v(0) 



= K*)7ov>7i> 
However, tak- 



while corrections are proportional to r\ 
ing into account these corrections is necessary to obtain 
a non-zero right-hand side of Eq. (TT2"]) . In finding these 
small corrections, we can employ the approximation of 
a "smooth distribution of local concentrations" used to 
proceed from Eq. ([28]) to (|30[) . that is, we can suppose 
the Tij values for all bonds ij of the given site i to be 
close to each other: 

Vij — fii — Tjj' (48) 

Then the solution of Eq. (|4T|) with the first-order correc- 
tions is provided by the relation: 



(49) 



Substituting into we again obtain Eq. 05J but 
correlator bfj in the effective mobility (f4T>]) is replaced by 



b h : 



M ij (b)=7^6?,expL5(A i + A i )/2]. 



(50) 



and 7®^ is the same as in (|41l) . 



Physically, the opportunity to reduce the vacancy- 
mediated kinetics to the equivalent direct exchange ki- 
netics is connected with the above-mentioned fact that 
in the course of evolution of an alloy, the distribution of 
vacancies adiabatically fast follows that of the main com- 
ponents. Therefore, one may suppose that such equiva- 
lence holds not only for simplified models (|3"6l or (|4"3"j) . 
but is actually a general feature of the vacancy-mediated 
kinetics, while for more general models, correlators bij, 
bfj or b% in Eqs. ([39]), (gBJ) or ([50]) are probably re- 
placed by some more complex expressions with similar 
properties. 

Function v(t) in Eq. ([41]) determines the rescaling of 
time between the initial VME model and the equivalent 
DAE model ([4"5"]) . Temporal evolution of this DAE model 
is actually described by the "reduced time" t r related to 
the real time t by the following differential or integral 
relations: 

dU = lf h dt = j av ~/h v v(t)dt; t r = f jf h (t')dt';(5l) 

Jo 

t= I ^ rf h (t' r )dt' r (52) 
Jo 

where = (7 Q ^)" 1 has the meaning of the mean time 
of an atomic exchange a ^ h, while the variable t r has a 
meaning of an effective number of such atomic exchanges. 
This natural physical variable is used below in describing 
the SSA simulation results. 
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To find the "rescaling function" v(t) in Eqs. (l36)) -(l52l). 
one should, generally, compare the results of simulation 
of precipitation based on the DAE model (|45|) to those 
based on the initial VME model. BV made such compar- 
ison for some simplified model of spinodal decomposition, 
while below we estimate v{t) for several realistic models 
of Fe-Cu alloys using comparison to the KMCA results. 
Note that the problem of rescaling of time between the 
real physical time and the time units employed in the 
simulation method used, e. g., number of Monte Carlo 
steps in the KMCA, persists in all simulations of VME 
kinetics (see Sec. IV C below), and it strongly depends, 
in particular, on the boundary conditions for vacancies 
adopted in simulations. For example, BV used the "va- 
cancy conservation" model, while in simulations [HI, H3] 
and below, a possible creation of vacancies at various lat- 
tice defects (grain boundaries, dislocations, etc) is taken 
into account. Thus, the form of the function v(i) de- 
pends also on the kinetic model used for vacancies. 

The results presented in Fig. [TOl below show that tem- 
poral variations of v{t) can be rather sharp. These vari- 
ations arise due to qualitative changes in the distribution 
of vacancies with respect to minority atoms related to the 
phenomenon of "vacancy trapping" at interfaces of pre- 
cipitates discussed in detail by BV and by Soisson and Fu 
17] (SF). The resulting excess of vacancies near growing 
or shrinking interfaces leads to an acceleration of effective 
exchange rates 7 efT with respect to the incubation stage 
when the precipitates are absent. It results in an increase 
of 7 eff (£) after be ginning of nucleation, and when the va- 
cancy trapping effect is strong, this increase can be very 
large, which is illustrated by Fig. [TOl below. At the same 
time, after the nucleation stage is over, the degree of this 
trapping does not change significantly. Therefore, the 
function v{t) can be expected to be approximately con- 
stant before and after nucleation and to monotonously 
increase with t in the course of nucleation, as illustrated 
by Fig. 3 of BV for v(t) in their simplified model. 



III. MAIN EQUATIONS OF STOCHASTIC 
STATISTICAL APPROACH 

A. Basic ideas of the classical theory of nucleation 

Before to describe the SSA, it is convenient to remind 
the main ideas of the classical theory of nucleation (CTN) 
[JQ. The CTN treats embryos of a new phase within the 
original metastable one as sufficiently large objects which 
arise due to thermodynamic fluctuations. The simplest 
version of the CTN considers the embryo as a homoge- 
neous droplet characterized by its radius R, the interface 
energy tr, and the free energy gain (with respect to the 
original metastable phase) per unit volume, Af. The 
excess free energy needed to form this embryo is: 

F{R) = 4irR 2 a - (4tt.R 3 /3) A/. (53) 



One of basic notions of the CTN is the critical embryo 
that can grow with no further loss of the free energy 
and thus with no fluctuations. For the model (|53p . it 
corresponds to the maximum of F(R) with respect to 
R, thus the critical radius R c and the nucleation barrier 
F c = F{R C ) are: 

R c = 2a /Af, F c = 16vro- 3 /3A/ 2 , (54) 

while the probability of the critical fluctuation needed to 
create this embryo is estimated according to the thermo- 
dynamic fluctuation theory 

W c ~ exp(-F c /T) - exp(-const a 3 /TAf 2 ). (55) 

Cahn and Hilliard [27| use d the Ginzburg-Landau-type 
free energy functional to allow for the diffuse character 
of the interface of the critical embryo, but their approach 
is valid only at high T ~ T c and for large embryos when 
the discrete lattice effects can be neglected. Dobretsov 
and Vaks [28|, [2t| developed a quantitative approach to 
calculate thermodynamics of critical embryos which takes 
into account the discrete lattice effects and uses the PC A 
rather than the simple MFA. Some results of this ap- 
proach are used below in Table II and Fig. O 

For supercritical embryos with R > R c , the CTN 
suggests fluctuation effects to be insignificant. There- 
fore, after completion of nucleation, the main type of mi- 
crostructural evolution is growth of embryos due to the 
diffusional flux of minority atoms from the matrix. Later 
on, the evaporation-condensation (or Lifshits-Slyozov- 
Wagner - LSW) mechanism becomes dominant when the 
larger precipitates grow at the expense of smaller ones 
0- Therefore, according to the CTN, decomposition 
of metastable solid solutions should include four well- 
defined stages (i) the incubation stage that precedes 
formation of first critical and supercritical embryos; (ii) 
the nucleation stage during which the supercritical pre- 
cipitate density d s reaches its maximum value; (iii) the 
growth stage when the density d s (t) remains approxi- 
mately constant but sizes of precipitates grow, and (iv) 
the coarsening stage when the density d s (t) decreases 
due to the LSW evaporation-condensation mechanism. 

These CTN ideas were confirmed by KMC simulations 
of Soisson and Martin [J] (SM) for some simple alloy 
model for which critical sizes and nucleation barriers (es- 
timated in Table II below) are rather large, see, e. g., Fig. 
1 of SM. However, in this work we mainly consider more 
realistic models of Fe-Cu alloys described in Sec. IV for 
which nucleation barriers and critical sizes are not large. 
Fluctuation effects in such alloys will be shown below to 
be strong and important not only for the nucleation but 
also for the growth stage. 

B. Stochastic kinetic equation and filtration of 
noise 

The quasi-equilibrium kinetic equations (1101) , (1381) and 
(1451) determine evolution of mean occupations of sites due 
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to the average atomic fluxes across each bond. However, 
these averaged equations can describe only those kinetic 
processes in which the total free energy F to t decreases 
(22l |23| , while the nucleation process should be accom- 
panied by a fluctuative increase of this F to t needed to 
overcome nucleation barriers. Therefore, to describe this 
process, we should consider fluctuations of atomic fluxes. 
The stochastic statistical approach for taking into ac- 
count such fluctuations was suggested by SPV (20|. In 
this section we present main equations of this approach 
while their refinements are described below in Sec. IV. 

We consider a binary alloy for which QKE has the form 
(l4"5|) . In the SSA, this QKE is replaced by the stochastic 
kinetic equation (SKE) which can be conveniently writ- 
ten in the finite difference form (for a short time interval 
St) given by Eq. (19) of SPV: 

Sc % = a{t + St) - Ci (t) = Scf + J2 Sn (y ( 56 ) 

Here Cj is the occupation of site i averaged over some 
locally equilibrated vicinity of this site, and the "diffu- 
sional" term Scf corresponds to the average atomic trans- 
fer to site i described by the right-hand side of the QKE 

Scf{c k } = My2sinh[/3(A 3 - A,)/2] St. (57) 
*»«(») 

The last term Sn{ 3 in the SKE ^ is the fluctuative 
atomic transfer through the bond (ij) which is described 
by the Langevin-noise-type method: each Sn{j is treated 
as a random quantity with the Gaussian probability dis- 
tribution: 

WiSnl) = A tJ exp[(-(5n/ i ) 2 /2A J -] (58) 

where Ay is the normalization constant, and the disper- 
sion Dij is the same as that for the actual fluctuative 
transfer Sn{j . This dispersion is related to the mobil- 
ity Mij and the time interval St in Eq. ([57} by the 
"fluctuation-dissipation" type relation (18) of SPV: 

D ij = ((Sn{ j ) 2 ) = 2M ij St, (59) 

Unlike standard applications of the Langevin-noise 
method to mechanical systems, for the non-uniform 
statistical systems under consideration Eqs. (f56|) - (j59|) 
should be supplemented by the "filtration of noise" pro- 
cedure that eliminates the short-wave contributions to 
fluctuations Snfj. As discussed in detail by SPV, these 
contributions to Eq. (|45p have been already included in 
the diffusional term Scf which is obtained by statisti- 
cal averaging just over these short-wave fluctuations. It 
agrees with the fact that all quantities entering Eq. ([56| , 
including the mean site occupation Cj, site chemical po- 
tentials \i, and the diffusional term Scf, have a physi- 
cal meaning only within some locally equilibrated region 



called in textbooks "a quasi-closed subsystem" [l| that 
contains a sufficiently large number of atoms. In other 
words, in our statistical description implying division of 
the whole non-uniform non-equilibrium alloy into locally 
equilibrated quasi-closed subsystems, we consider only 
"thermodynamic" fluctuations Sn{j which have approxi- 
mately the same value (for all bonds ij of the given crys- 
tal orientation a) within each quasi-closed subsystem, 
while a non-zero fluctuative contribution to the total Sci 
in (|56p arises only due to a relatively weak non-uniformity 
of these fluctuations. Therefore, in the last term of Eq. 
(f56")) . the full fluctuative transfer Sn{^ should be replaced 

by its long-wave (or "coarse-grained") part Sn{^. The 
latter can be obtained by introducing a proper cut-off 
factor F c (k) in the Fourier-component <5n/ Q (k) of the 
full fluctuation 6n{j = (5n^(R SQ ) where R SQ denotes the 
position of the bond ij center in the appropriate crystal 
sublattice a formed by these centers [201 ] : 

Sn{ c (H sa ) = ^exp(-ikR SQ ) Sn f a (k) F c (k) 

k 

Sn f a (k) = jrYl ex P(* kR -) 8n f a (R sa ) (60) 

where N is the total number of lattice sites (or atoms) 
in the crystal. The cut-off factor F c (k) can be taken in 
the simple Gaussian-like form which for the BCC lattice 
looks as follows: 

F c BCC (k) = exp [-4g 2 (l - cost^i cos ip 2 cos tp 3 )] (61) 

where ip v — k v a/2; k v is the vector k component along 
the main crystal axis v\ and a is the BCC lattice con- 
stant. 

At large g 2 » 1, the expression (|6T|) is reduced to 
a Gaussian exp (— k 2 l 2 /2) with / = ga. Thus, the re- 
duced length g = l/a characterizes the mean size of lo- 
cally equilibrated quasi-closed subsystems. Generally, it 
is the characteristic length of uniformity of site chemi- 
cal potentials Ai, which for a single-phase state (before 
nucleation) coincides with the uniformity length for local 
concentrations Estimates of this size are discussed 
below in Sec. IV B. 



IV. MODELS AND METHODS OF 
SIMULATIONS 

A. Alloy models and states used for simulations 

In the most of our simulations we use the first-principle 
microscopic model of Fe-Cu alloys suggested by SF [17j 
that describes well available thermodynamic and kinetic 
data for such alloys. This is the second-neighbor interac- 
tion model based on ab initio Density- Functional Theory 
calculations using SIESTA code. For this model, config- 
urational interactions v n for the n-th neighbors in (|22l) , 
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FIG. 1: (color online). Calculated phase diagrams for the SF 
model (upper figure) and LBS moel (lower figure) described in 
the text. Thick lines (red online) correspond to the pair clus- 
ter approximation (PC A), and thin lines (green online) cor- 
respond to the mean-field approximation (MFA); solid lines 
are binodals, and dashed lines are spinodals. Solid circles and 
chained line show the binodal and spinodal calculated using 
the CALPHAD expression for free energy of bcc Fe-Cu alloys 
taken from [ll[. Triangles indicate the (c, T) states used for 
simulations. 



kinetic interactions = u n in (|17p . saddle-point en- 
ergy parameters A p in ([55]) . activation energies E^ in 
([TJJ (all in eV) and attempt frequencies w pv in (fT4")) (in 
sec" 1 ) have the following values: 



vi = -0.121 + 0.182T, 
ui = 0.127 - 0.091T, 
A Cu = 0.05, A Fc = 
E^ v = 0.438, 



v 2 = -0.021 + 0.091T, 
it 2 = 0.044 - 0.045T, 
-0.03 



wcuv = 5-10 



15 



El™ = 0.698, 
ojpcv = 2-10 



15 



(62) 



Here the second terms in expressions for v n and u n de- 
scribe phenomenologically the influence of anharmonic 
and magnetic effects. Our simulations for the SF model 
were mainly performed at concentration c = 0.0134 and 
three temperatures T: 773, 713 and 663 K; these alloy 
states will be denoted as the SF-1, SF-2 and SF-3 state, 
respectively. In Tables I and II and in Sec. V we also 
present some results for the alloy state SF-4 that corre- 
sponds to c = 0.0197 and T=773 K. 

We also employ two more models suggested by Le 
Bouar and Soisson [l6[ (LBS). Parameters of these mod- 
els have been fitted to energies of configurations and va- 
cancy migration barriers computed with the Embedded 
Atom Method potential by Ludwig et al. [3(| ■ The model 
LBS-1 is the nearest-neighbor interaction model with the 
following parameter values (in the same units as in (|62l) ): 



LBS - 1 



vi = -0.20, ui = 0.155, 

A Cu = A Fc = q 



£ a c c uv = 0.018, 



ElT = 0-64, 



wcuv = w Fov = 5-10 



(63) 



For the model LBS-2, all parameters are the same as in 
(|6"3")) except for the saddle-point-energy parameter A Fc 



larger sizes of nucleating precipitates), and the larger dif- 
ference between vacancy formation energies in pure BCC 
iron and pure BCC copper (which leads to a stronger 
vacancy trapping in copper precipitates and to a higher 
mobility of precipitates). As discussed by SF, the LBS 
models describe Fe-Cu alloys less realistically than the 
SF model, but we will consider these LBS models for 
methodical reasons, to follow the influence of variations 
of interaction constants and saddle-point energy param- 
eters on the precipitation kinetics. 

Binodals and spinodals for the SF and LBS models 
calculated using the PCA , MFA and CALPHAD ex- 
pressions for free energy are presented in Fig. [T] where 
we also show the alloy states used in our simulations. 
These states are chosen in the metastable region T s (c) < 
T < Tb(c) or Cb(T) < c < c s (T) (where index s or 
b corresponds to the spinodal or binodal) which corre- 
sponds to the nucleation and growth (NG) type of alloy 
decomposition. The degree of supersaturation for each 
of these states can be quantitatively characterized by the 
reduced supersaturation parameter s introduced in [2^ | : 



s(c,T) = [c-c b (T)]/[c s (T)-c b (T)]. 



(65) 



Values s < 1 correspond to the NG, while s > 1, to 
the spinodal decomposition (SD) evolution type; for the 
alloy states studied, values of s are presented in Table I. 

To appreciate accuracy of results presented in Fig. [1] 
we first note that the PCA fully takes into account the 
pairwise correlations of atomic positions and neglects 
only many-particle ones [26|. Therefore, in the dilute 
alloy limit c$ — > 0, the PCA expressions for thermody- 
namic potentials become exact. In particular, the free 
energy (|24[) in this limit takes the form 



"PCA 



c->0 



2 ^ 



CiCj fij j 



(66) 



which is just the discrete lattice analogue of the first term 
in the virial expansion of the free energy in powers of 
density (in our case, in powers of Ci) for the standard 
"weakly non-ideal gas" theory [l|. Thus, at low concen- 
trations under consideration, we can expect the accuracy 
of all thermodynamic results of the PCA, including both 
binodals and spinodals, to be high, and for all thermo- 
dynamic calculations we use the PCA. 

TABLE I. Alloy states used in our simulations. 



LBS - 2 : A Fc = -0.238 (64) 


Model 


SF 


LBS-1,2 


Our simulations for these two models were performed at 


T, K 


773 


713 


663 


773 


1000 


c=0.01 and T=1000 K, and these alloy states will be 














denoted as the LBS-1 and LBS-2 state, respectively. 


c, at% 


1.34 


1.34 


1.34 


1.97 


1 


Comparison of SF and LBS models discussed in [TtJ 














(a more detailed comparison is given in Ref.(3lJ) shows, 


s 


0.285 


0.352 


0.425 


0.426 


0.459 


that the SF model yields a lower solution energy of cop- 
per in BCC iron (which leads to the lower density and 


Alloy state 


SF-1 


SF-2 


SF-3 


SF-4 


LBS-1,2 
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FIG. 2: (color online). Concentration profiles Ad = c; — c 
in the critical embryos for the alloy states studied. Different 
curves from top to bottom correspond to the state SM 
SF-1, SF-2, SF-4, SF-3 and LBS, respectively. 



Fig. Q] shows that at relatively low c and T values 
used in our simulations, the binodals calculated using 
the PCA, MFA or CALPHAD methods do virtually co- 
incide with each other. At the same time, the CAL- 
PHAD or MFA-calculated spinodals c s (T) exceed those 
found in the PCA (and thus, the exact ones) by about 
twice. It leads to a strong underestimating of the re- 
duced supersaturation s(c, T), which, in its turn, should 
result in drastic distortions of microstructural evolution. 
Therefore, as mentioned in Sec. I, using the phase-field 
type methods based on CALPHAD or MFA expressions 
for thermodynamic potentials can hardly provide an ad- 
equate description of alloy decomposition kinetics at low 
concentrations and temperatures under consideration. 

Fig. [2] illustrates the structure of "thermodynamic" 
critical embryos for the alloy states studied. This struc- 
ture was calculated by the method of Ref . [28| mentioned 
in Sec. Ill A with the use of the PCA. The figure shows, 
in particular, that the sizes of critical embryos in our 
problem are comparable to the host (BCC iron) lattice 
constant: ap c = 0.287 nm. In Table II we present main 
characteristics of these critical embryos: the reduced nu- 
cleation barrier F c /T, the mean radius of the embryo, 
i? c , the total excess of minority (copper) atoms in the 
embryo with respect to the initial state, AN C , and the 
total number of minority atoms in the embryo, N c . The 
barrier F c is calculated according to Eq. (4) in (29j 
(where it was denoted as AQ C ), while quantities R Cl 
AN C , and N c are defined by relations: 



TABLE II. Parameters of "thermodynamic" critical 
embryos for the alloy states considered. 



All j- j- 

Alloy state 


s 


F c /T 


i? c , nm 


AN C 


N c 


niv t r A 1 

SM |4J 


0.287 


7.48 


0.488 


32.5 


33.8 


SF-1 


0.285 


4.38 


0.417 


13.3 


14.3 


SF-2 


0.352 


2.47 


0.438 


9.7 


10.7 


SF-3 


0.425 


1.36 


0.468 


7.2 


8.2 


SF-4 


0.426 


1.90 


0.475 


10.4 


11.9 


LBS-1,2 


0.459 


0.83 


0.497 


5.4 


6.2 



Let us now discuss the expressions for effective mo- 
bilities Mij in Eq. (|4"5j) . Model LBS-1 corresponds to 
the configuration-independent saddle-point energies for 
which correlators in (11 Tp do not depend on the kind 
p of an atom. As mentioned in Sec. II C, for this case 
the adiabaticity equation (|34f is solved analytically, and 
the effective mobility M l3 = M% h is given by Eq. (|39]) . 

For two other models, SF and LBS-2, the saddle-point 
energies depend on configurations. Thus the adiabaticity 
equation (|34[) can be approximately solved only if either 
of inequalities (|43j) is obeyed at Cj values significant for 
the kinetic process studied; for brevity we denote such Cj 
as c s . Let us first consider NG-type processes for the SF 
model. The experience of our simulations for all models 
considered (illustrated by Figs. [20ll22l below) shows that 
the "significant for NG" local concentrations c s have 
usually the order c s ~0.1. To estimate products rjiTij in 
Eqs. ([4"3"]). we first consider the case when these c s are 
small, thus both the functions ln(l— QijCj) ~ ln(l — fijC s ) 
inEq. flUJ) and (1 + c%/ Ap ) 6 - exp[6 ln(l + c s / Ap )] in 
Eq. (|3"3"|) can be expanded in powers of c s . Substitut- 
ing numerical values of parameters (|62p into these expan- 
sions we obtain the following estimates of £y = r\itij for 
the SF-1 and SF-3 states: 



i i 

N c = AN C + cN^. (67) 



Here index i numbers the lattice sites, N£ is the total 
number of sites (atoms) in a precipitate, and the term 
cN% describes the contribution to N c of the "uniform 
background" minority atoms. For comparison, in Fig. 3 
and Table II we present also the analogous results for 
the SM Q alloy state which corresponds to the nearest- 
neighbor interaction v\ < 0, T — 0.4(— vi), and c = 
0.03. 



4 s /" 1 - 20c s exp (80c s ), ^f/" 3 ~ 40c s exp (160c.) 

(68) 

while for the SF-2 state, the £,j value lies between these 
two estimates. For c s >0.1, all these values much exceed 
unity. Therefore, inequality (b) in (|4"3")) is obeyed, and 
the QKE (|45|) with the effective mobility My given by 
Eq. ([5T?]) can be used for all our SF states. 

To estimate products £y = TjiVij in Eqs. (I43|) for the 
LBS-2 state, we again expand functions ln(l — gijCj) ~ 
ln(l - f ijCs ) in ||2DJ) and (1 + c7j/ Ap ) 6 - exp [6 In (1 + 
c s / Ap )] in (|3"3"|) in powers of c s . Substituting numerical 
values of parameters from Eqs. (|63[) . we obtain in this 
case: 

?7z^ BS ~ 2 - 1400 c s exp (-40c,). (69) 
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At c s between 0.01 and 0.1, the right-hand side (rhs) of 
this estimate is large, inequality (b) in (|4"3")) is obeyed, 
and Eq. (|50|) for the effective mobility My can be used; 
at c s >0,15, this rhs is small, and Eq. (|4"6")) for can 
be used; and at c s between 0.1 and 0.15, the rhs is of the 
order of unity, thus neither of inequalities (|43| is obeyed. 
Therefore, the equivalence of the VME kinetics to that 
for the DAE model (|45p can not be formally proved for 
the LBS-2 state, and we made no the SSA simulations for 
this state. However, the above-mentioned remarks enable 
us to suggest that actually such equivalence probably ex- 
ists, with the mobility in the QKE (|4"5)) smoothly 
varying between expression (|50[) for c s <0.1 and expres- 
sion (|46|) for c s >0.15. Then the precipitation kinetics 
for the LBS-2 state should be quite similar to that for 
the LBS-1 state differing only by replacing in Eq. (|4"6"]) 
the factor bfj = bij for the LBS-1 state by some smaller 
factor varying between &^<6,j at c s < 0.1 and bfj = fey 
at c s >0.15. It should correspond just to some slowing 
down of precipitation kinetics for the LBS-2 state with 



respect to the LBS-1 state, and the KMCA results of 
LBS [l|| seem to agree with these considerations. 

To characterize strength of configurational (or "ther- 
modynamic" ) interactions v n and kinetic interactions 
u n for the alloy states considered, in Table III we present 
values of Mayer functions and of analogous "kinetic" 
functions f% defined by relations: 

II = exp (-(iv n ) - 1, ft = exp (/3u n ) - 1, (70) 

as well as functions /£ in Eq. ([55]) . Functions enter 
Eqs. (l20l) - (|24p for site chemical potentials Ai and the free 
energy F, while functions and ft enter Eqs. (|25p- 
(1551) . (|3"9"j) . pp| andEHl) for correlators b\ } and effective 
mobilities My. In Table III we also present expressions 
for correlators bfj' and for reduced effective mobilities 
M(j = Mijhf h in Eqs. (®, dMD and at small 
Cj ~ Cij < which have been mentioned to be most 

significant for the NG-type processes. 



TABLE III. Values of functions f%, f%, / A in Eqs. ([701) and (1551) and expressions for correlators bfj and reduced 
effective mobilities M[j = My/7^ in Eqs. ([55]) . (|4l)|) and (|5D|) at small q ~ cy < for the alloy states studied. 



Alloy state 


fv 

Jl 


£V 

11 


fu 

Jl 


fu 
J2 


fFc 

J a 


fCu 
J A 


6 Fc 
«.? 


Mr. 

« 


SF-1 


4.1 


0.3 


5.2 


0.9 


5.8 


24 


exp (81^) 


exp (47c») 


SF-2 


4.9 


0.3 


6.3 


1.0 


7.0 


33 


exp (98q) 


exp (57cj) 


SF-3 


5.9 


0.3 


7.5 


1.1 


8.5 


43 


exp (118cj) 


exp (69cj) 


LBS-1 


9.2 





5.0 





5.0 


5.0 


exp (71c,) 


exp (-3c,) 



Let us discuss the results presented in Table III. First, 
they show that both the thermodynamic and kinetic in- 
teractions for the alloy systems considered are rather 
strong: the /J, /" and / A values much exceed unity. 
It again shows that the MFA or CALPHAD-type expres- 
sions for thermodynamic and kinetic parameters based 
on the approximations (3v n 1, (3u n -C 1 [2l| can not 
be used to describe these alloy states. Second, the last 
column of Table III shows that the NG kinetics for the SF 
model should notably differ from that for the less realistic 
LBS-1 model. At small local concentrations c, consid- 
ered, the reduced mobility M[j for the LBS-1 model is 
virtually a constant close to unity, while for the SF model 
it sharply rises with Cj and is typically very large. Ac- 
cording to Eqs. ([58]) and ([59]) . this mobility determines 
scale of fluctuativc terms in the SKE ([56]) . Therefore, we 
can expect the manifestations of fluctuation effects for 
the SF model to be much stronger than for the LBS-1 



model. It agrees with the KMCA results presented be- 
low in Figs. Eni 



B. Estimating the local equilibrium length for the 

SSA 

As discussed in Sec. Ill B, the reduced length I = ga 
in the SSA equations (|60| and (|6TT) characterizes sizes of 
locally equilibrium quasi-closed subsystems used in our 
statistical description of a nonequilibrium alloy. This 
length can not be chosen lower than the characteristic 
length of non-uniformity of local chemical potentials, ^ nu , 
which for the nucleation processes has typically the same 
order of magnitude as the critical embryo size R c , see, 
e. g., Figs. I201I221 below. The actual distribution of local 
equilibrium lengths I < l nu in an alloy varies with both 
space and time; in particular, after creation of a super- 
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FIG. 3: Total number of precipitates that contain i > p cop- 
per atoms, N p (t r ,g), versus the reduced time t r defined by 
Eq. (|51[) obtained in SSA simulations with different g for the 
SF-1 state. Frames (a), (b) and (c) correspond to g equal to 
1.75, 1.65 and 1.55, while different curves from left to right 
in each frame correspond to p = 10, 20, 40, 60 and 70, respec- 
tively. 

critical precipitate, the degree of local equilibrium in the 
adjacent region should significantly increase with respect 
to other regions where such precipitates are not born yet. 

For simplicity we will characterize the distribution of 
all local lengths I by a single spatially averaged parame- 
ter I = ga where the reduced length g, generally, varies 
with the evolution time t or the reduced time t r (|5ip 
used in the SSA. After completion of nucleation at some 
t r = , the alloy rapidly approaches the full two-phase 
equilibrium, and the length / should become large. Then 
the cut-off parameter g = g(t r ) in Eq. (f6"Tj) at t r > t^f 
should be large, too, fluctuative terms Srif in (|55| be- 
come small, and the SKE (151)1) transforms into the QKE 
(|45|) with no fluctuation terms. 

To describe the above-discussed physical picture with 
the minimal number of model parameters, we approxi- 
mated the time dependence g(t r ) by the following sim- 
plest expression: 

t r < : g(t r ) = g a 

tr>t?: g 2 (t r )=g 2 + (t r -t?)CD cS . (71) 

Here D is the effective reduced diffusivity which for the 
direct-exchange model (|45)l can be estimated as D eS ~ 
Tf^Cu dHHH; and C is a numerical factor, e. g., C =2 
for the standard diffusion law. At C ~ 1, the evolution 
of microstructure was found to virtually not depend on 
the C value, and usually we put C =2. 

The simple model (jTTj) for g(t r ) includes an unphysical 
break at t r — t^f which leads to the presence of analo- 
gous fictitious breaks in various characteristics of evolu- 
tion, e. g., in Figs. flUTI 151 fTUl and [TBI As mentioned, in 
reality the effective equilibrium length starts to increase 
immediately after beginning of nucleation, thus the func- 
tion g(t r ) — l/a monotonously increases with time with 
no breaks. Therefore, for a more realistic description, we 
should use for g(t r ) some other models which describe 
its continuous increase with t r and the resulting smooth 
decrease of fluctuations in the course of both nucleation 
and growth stages. However, the experience of our simu- 
lations shows that at t r not close to t^f , the main char- 
acteristics of microstructure, such as the density and sizes 
of supercritical precipitates, are not sensitive to the de- 
tailed form of g(t r ) provided the scale of this function 
is determined by the "maximum thermodynamic gain" 
principle described below. Thus, in this work we employ 
for g(t r ) the simplest form (f7T|) . 

The parameter go in Eq. (|71[) can be estimated by 
two ways. First, it can be found by fitting the SSA sim- 
ulation results for the evolution of density of precipitates 



to the analogous KMCA results, for example, to those 
presented in Figs. ITlTHl below. However, such "KMCA- 
based" estimates of go would restrict possible applica- 
tions of the SSA by the models for which reliable KMCA 
results are available. To estimate go within the SSA, 
we can try to extend the second law of thermodynamics, 
that is, the principle of the free energy minimum with re- 
spect to all its free parameters valid for equilibrium sys- 
tems, to the kinetic processes in non-equilibrium systems 
studied. To this end we note that the main characteris- 
tics of microstructure formed in the course of the nucle- 
ation process, such as the characteristic non-uniformity 
length l nu ~ I mentioned above, can be considered as 
"free" parameters of a nonequilibrium state analogous to 
"static" free parameters for equilibrium systems. Thus 
it seems natural to suggest that the kinetic path of evo- 
lution of this nonequilibrium state should correspond to 
the maximum thermodynamic gain, that is, to the max- 
imum rate of decrease of free energy. This suggestion 
can extend the "excess entropy production" approach 
to thermodynamics of irreversible processes discussed by 
Prigogine and coworkers [34| to the kinetics of essentially 
non-uniform and non-equilibrium systems under consid- 
eration. Then the characteristic value l nu ~ god can be 
estimated from the condition of the maximum thermody- 
namic gain in the course of the nucleation process, which 
in our model (|7ip corresponds to the minimum of the 
free energy F(g ,t r < t^f) with respect to go- 

This maximum thermodynamic gain should evidently 
correspond to the formation of maximum number of large 
supercritical precipitates. At the same time, the density 
dp of such precipitates sharply depends on the effective 
size I of quasi-equilibrium systems for the nucleation pro- 
cess. When g = l/a is too large, the fluctuations are too 
weak to overcome nucleation barriers, while when g is 
too small, the fluctuations are too strong to allow forma- 
tion of large and steadily growing precipitates. There- 
fore, the dependences d p (g) should have a pronounced 
maximum at some optimal g = go- These considerations 
are illustrated by Fig. [3] where we present temporal de- 
pendences of the total number of precipitates containing 
i > p copper atoms, N p (t r ,g), obtained in the SSA sim- 
ulations with different g. The "end of nucleation" time 
t~? here and below is defined as the time of creation of 
the last "critical" precipitate with p > p c , while these 
p c values (which for SF alloy states are close to the N c 
values in Table II) are estimated in Sec. V A. Such defini- 
tion of is somewhat arbitrary but it makes virtually 
no effect on the simulation results. For the evolution 
shown in frame [3^, the fluctuations seem to be too weak, 
thus the length I — ga is too large. For frame [3};, the 
fluctuations are evidently too strong which leads to an 
unphysical evolution of large precipitates: they do not 
grow, in disagreement with the second law of thermody- 
namics; thus, the g value here is too small. Finally, for 
the evolution shown in frame (3Jd, the g value seems to 
be close to optimum. 

In Figs. HE we present temporal dependences of the 
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FIG. 4: (color online). Temporal dependence of the free en- 
ergy per copper atom, F(t r ,go), obtained in SSA simulations 
with different go for the SF-1 state. Curves 1, 2, 3. 4 and 
5 (red, green, blue, purple and black online) correspond to 
go =1.55, 1.65, 1.7, 1.75 and 1.85, respectively. 

FIG. 5: (color online). Same as in Fig. 2] but for the SF- 
2 state. Curves 1, 2, 3 and 4 (red, green, blue and purple 
online) correspond to go =1.6, 1.7, 1.8 and 1.9, respectively. 

free energy per copper atom, F(t r ,go), obtained in the 
SSA simulations with different g in Eq. flTTj) . For 
deiniteness, the initial state for these simulations was 
taken uniform: Ci(0) = c=const, thus the initial increase 
of F at t r < 0.1 1^: is related just to switching-on fluc- 
tuations at t r = and has no physical meaning. At 
t r ~ t^f , functions F(t r ,go) show the above-mentioned 
fictitious breaks due the analogous breaks in our model 
function g(t r ) in (|7"Tj) . However, at t r not close to zero 
or to t^, say at t r ~ 0.5^, functions F(t r ,go), pre- 
sented in Figs. 2][7] can be realistic. Therefore, the "opti- 
mal" go value can be estimated from comparison of these 
functions at t r ~ 0.5^ for different go . 

Figs. 0][7]show that these functions have a distinct min- 
imum at certain go for each alloy state considered. For 
the SF-1 state, this minimum corresponds to go equal 
to 1.65 or 1.7; for the SF-2 state, to g —1.7 or 1.8; 
for the SF-3 and LBS-1 states, to g Q = 1.9 or 1.8. The 
first values seem to be a bit more appropriate, but us- 
ing in simulations the second values changes results only 
slightly; it is illustrated below for the SF-2 model. There- 
fore, for the reduced length go in Eq. (fTTj) we use the 
following values: 

fig*" 1 = 1.65; 5o SF - 2 = 1.7; 

ffo SF - 3 = 1.9; flj 188 " 1 = 1.9. (72) 

The minimum local equilibrium lengths lo — god for 
these go have the same order of magnitude as critical 
sizes R c in Table II, in accordance with the considera- 
tions mentioned above. 

Let us comment on the loss of validity of the SSA at 
low g which is manifested, in particular, in the above- 
mentioned unphysical results presented in frame This 
problem was discussed in detail by SPV [20j who noted 
that the statistical approach used, in particular, ba- 
sic equations (l9"))-([l"3"f that include averaging over lo- 
cally equilibrated quasi-closed subsystems, imply their 
reduced size g to be not too small, so that they include 
a sufficiently large number of atoms, while site chemical 
potentials Ai within these subsystems should obey the 



FIG. 6: (color online). Same as in Fig. [4] but for the SF- 
3 state. Curves 1, 2, 3 and 4 (red, green, blue and purple 
online) correspond to go =1.8, 1.9, 2 and 2.1, respectively. 



FIG. 7: (color online). Same as in Fig. [I] but for the LBS- 
1 state. Curves 1, 2, 3. 4 and 5 (red, green, purple, blue 
and black online) correspond to go =1.7, 1.8, 1.9, 2 and 2.1, 
respectively. 

FIG. 8: (color online). Parameter of non-equilibrium (|73[) for 
the incubation stage observed in SSA simulations with differ- 
ent reduced lengths g. Circles and squares (red and black 
online) correspond to the SF-1 and LBS-1 state, respectively. 
Arrows indicate go values in Eqs. (|72l) . 

local equilibrium condition Ai ~ const. The scale of vio- 
lations of this basic condition can be characterized by the 
"parameter of non-equilibrium" J introduced by SPV: 

where Nf, is the total number of the nearest-neighbor 
bonds {ij}, and the sum is taken over all such bonds in 
an alloy. In Fig. [5] we show the values of this param- 
eter averaged over incubation stage, Ji nG (<?) = (J)incl 
for different alloy states, functions Ji nc (g) are similar. 
At small g < 1.5, these functions start to sharply rise 
which reflects sharp violations of statistical equilibrium 
within too small quasi-closed subsystems. However, for 
g=go values presented in Eqs. ([T2"|). this parameter is 
still small: Jmc(ff) ~ 0,1-0.15, thus employing the SSA 
seems to be justified. 

In Fig. [9] we show temporal dependences of the 
nonequilibrium parameter (I73[) at first stages of evolu- 
tion. As mentioned, breaks in curves J (t^ ) at t r — t r 
are due to the similar breaks in our simple model (|71[) . 
while for more realistic models mentioned above these 
breaks should be replaced by some smooth decrease of 
J(t r ) at t r >0.5t^. The increase of J(t r ) after begin- 
ning of nucleation is related to arising of interfaces and 
spacial non-uniformities which lead to some additional, 
"non-uniform" contributions to differences | A^ — Xj | in 
Eq. 03]). 

C. Methods of KMCA simulations 

The KMCA used in this study is described in detail in 
Refs. [H, [I3|- We just recall here the physical principles 
of the method, underlining the difference with the SSA. 
The KMC simulations follow the evolution of the atomic 
configuration in a simulation box of N = 64 3 lattice sites 
containing Fe and Cu atoms and one vacancy, with peri- 
odic boundary conditions. At each Monte Carlo step, 8 



FIG. 9: (color online). Parameter J(g,t r ) (I73[) observed in 
the SSA simulations with g—go from (|72l) . Curves 1, 2, 3 
and 4 (red, green, blue and purple online) correspond to the 
SF-1, SF-2, SF-3 and LBS-1 alloy state, respectively. 
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atom-vacancy exchanges between nearest-neighbor sites 
can occur in a BCC lattice, with the jump frequencies 
W[ 3 V given by Eq. ((6|). The activation energies are 
exactly computed for each local configuration, without 
using any mean-held approximation. One of these ex- 
changes is chosen, using a pseudo-random generator, by 
means of a residence time algorithm (|35{. The physical 
time of the Monte Carlo simulation is given by: 



tuCS — 1 , 



(74) 



where the sum runs over the 8 possible jumps. This time 
must be rescaled to take into account the real vacancy 
concentration, which depends on the precipitation mi- 
crostructure. If one assumes that the vacancy concentra- 
tion remains at equilibrium in the different phases during 
all the precipitation, a convenient way to perform this 
time rescaling is : 



„MC 

t = tMCS—eq- 



(Fe) 



(Fe) 



(75) 



where Cy (Fe) is the equilibrium vacancy concentration in 
pure iron and Cy (Fe) is the vacancy concentration in 
the copper-free regions of the KMC simulation box [lH 
Il7| . The number N p of copper-rich precipitates and their 
average size (i) (discussed in Sec. V A) are computed 
by considering only clusters which contain i > p copper 
atoms connected by at least one nearest-neighbor bond. 



D. Methods of SSA simulations 

All SSA simulations were performed in the cubic box 
containing N=2 x (64) 3 BCC lattice sites with periodic 
boundary conditions. In describing precipitates, the pre- 
cipitate containing i > p copper atoms is defined as a set 
of adjacent sites i (connected by at least one bond) for 
which their mean occupations Cj exceed a certain cut-off 
value: c, > c cu t, while this set contains not less than p 
copper atoms: 



N, 



Cu 



z2 Ci > p- 



(76) 



The choice of c cu t was found to be not essential, and 
for definiteness we took it the same as in experimental 
studies 0: c cut =0.05. 

In solving the SKE (|5"oT) with the diffusional term (1571) , 
we should take into account that this term is propor- 
tional to product of the generalized mobility My and 
the factor 2 sinh[/3(Ai — Xj)/2] ~ /3(Aj — Xj) describing a 
thermodynamic driving force, while mobilities My given 
by Eqs. J3H|) or ([ST)]) are proportional to the correlator fog- 
or by very sharply rising with the local concentrations 
Cjj the latter is illustrated by Eq. and by two last 
columns of table III. These very sharp dependences do 



not allow us to use for solving the SKE the standard it- 
erative methods, such as the Euler or Runge-Kutta ones: 
after several iterations, the product My/3(Aj — Xj) be- 
comes so large that the time step needed to achieve a 
numerical stability gets too small for these algorithms 
can be used. On the other hand, there is no physical rea- 
son for the diffusional term to be too large, as the high 
mobility My should lead to a very fast approaching the 
local equilibrium state at which local chemical potentials 
of adjacent sites, Xi and Xj, are very close to each other. 
Thus, in reality diffusional terms (|57p remain to be rea- 
sonably small. The problem with application of standard 
iterative methods arises due to their discrete nature, thus 
numerical values of differences |A,- — Xj\ can not catch up 
with a very fast increase of 6y , and the fictitious increase 
of their product happens. 

To overcome this methodical difficulty, in our itera- 
tive computations we put restrictions on fry setting it 
to not exceed a certain value b"* ax : < b™; ax . Then 
we made simulations with different b"] ax , from smaller 
to larger ones, until all physically important characteris- 
tics of evolution, including the incubation and nucleation 
time, i*™ c and t^f, and the maximum density of super- 
critical precipitates, d™ ax , ceased to significantly change 
under increase of b™ ax . It happens at b™- ax = 250, and 
this value was used in all our simulations. Actually, at 
6y ax — 500, the d™ ax values may be even lower than at 
b rnax = 25Q (while for b^ ax < 250, they monotonously 
increase with b™; ax ) but differences lie within statistical 
errors. After the nucleation stage is over and fluctua- 
tions are switched off according to model (|7Tj) , values 
|/3(A<j — Xj | typically decrease by two orders of magnitude 
with respect to the nucleation stage. This allows us to 
significantly decrease the b™- ax usually up to ft-™ 02 = 10, 
and thus to increase the time-step in solving equations. 
Again we checked that using the same b™j ax equal to 250 
both before and after t^f , and reducing it at t r > 



up to W 



=10, lead to the same description of evolu- 



tion. In our case, using the Runge-Kutta method does 
not provide any improvement of stability for the numeri- 
cal solution of equations while it needs more calculations 
at each integration step, so in our simulations we used 
the more simple Euler method. 

Finally, let us discuss the rescaling of the reduced time 
t r used in the SSA to the physical time t determined by 
Eq. (|52[) . In accordance with the remarks in the end of 
Sec. II C, we suppose the effective mean time of direct 
atomic exchanges, in (fB"2"]) . to be constant before and 
after nucleation (that is, at both t r < t) nc and t r > t^f), 
and to linearly vary with t r in the course of nucleation: 



bf \ Lj, , 



r ah = < a 1 + (a 2 - ai ) / N 'J. , C c < U < t?; 



(t? - t r nc ) 



«2. 



t r > t 



N 



(77) 

The constant a\ is the ratio of physical and reduced in- 
cubation times: a± = Unc/fJ 10 , while the constant a 2 
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FIG. 10: (color online). Temporal dependence of effective 
direct exchange rate 7cuf c in Eq. (|4ip estimated from com- 
parison to the KMC A results as described in the text. Differ- 
ent curves from top to bottom correspond to the alloy states 
LBS-1, SF-1, SF-2 and SF-3, respectively. Arrows indicate 
the incubation time ij ne . 



FIG. 11: (color online). Evolution of density of precipitates 
containing i > p copper atoms, d p (left scale), and the to- 
tal number of such precipitates within the KMC simulation 
box V S KMC , N p = N* MC (right scale), for the SF-1 alloy 
state. For the SSA simulations with 1/ S SSA =2V S KMC , these 
N p should be doubled: Np SA =2N p . Solid lines from top 
to bottom (green, red, blue and pirple online) correspond to 
the KMCA results for p = ll, 15, 21 and 26, respectively, 
while symbols (circles, squares, triangles and crosses) show 
the analogous results obtained in another KMC run. Dashed 
line shows the SSA results for p = 15, go = 1.65. 

can be estimated from the fit of the SSA simulation re- 
sults to the evolution rate at the coarsening stage. Thus, 
model (|77p includes only two parameters, U nc and a^, 
which can be estimated either from comparison to KMC 
simulations or from experiments. 

In this work we use for such estimates the KMCA re- 
sults described below. The resulting rescaling of time is 
presented in Fig. [TO] as the dependence of effective direct 
exchange rates 7°^ in Eq. (|5ip (to be abbreviated 7 cff ) 
on the physical time t; this dependence has a more clear 
physical meaning than T a h{t r ) in Eq. (1771) . The results 
presented in Fig. [TO] and in analogous Fig. 3 of BV [22| 
clearly illustrate the decisive role of vacancy trapping for 
temporal dependences of effective direct exchange rates. 
For all models considered, these rates monotonously in- 
crease with the evolution time t which reflects the de- 
velopment of vacancy trapping in the course of precipita- 
tion. For the BV model, this trapping is relatively weak, 
thus 7 cff increases to the coarsening stage just by about 
twice. For the LBS model, the vacancy trapping is also 
not too pronounced (though stronger than for the BV 
model), thus the full increase of 7 cff is about 6 times. 
For the more realistic SF model, the vacancy trapping is 
very strong [l7[ . Thus for all three SF states considered, 
the 7 cff values rise between the incubation and nucle- 
ation stages by more than two orders of magnitude, and 
with lowering temperatures T this rise seems to increase, 
in accordance with a probable more strong trapping at 
lower T. 



V. EVOLUTION OF MICROSTRUCTURE 
OBSERVED IN KMCA AND SSA SIMULATIONS 

A. Evolution of density and sizes of precipitates 

Evolution of density of different precipitates is shown 
in Figs. [TTllTil Solid lines and symbols in these figures 



FIG. 12: (color online). Same as in Fig. [TT]but for the SF-2 
state and f> = 8, 11, 16 and 21, respectively. Dashed curve 
shows the SSA results for p = 11, go — 1.7, and chained 
line, the SSA results for p — 11, go = 1.8. 



FIG. 13: (color online). Same as in Fig. [TT]but for the SF-3 
state and p = Q, 8, 11, 16 and 21, respectively, while dashed 
curve shows the SSA results for p — 8, go — 1.9. 

show the results obtained in two different KMC runs; 
differences between these lines and symbols illustrate the 
scale of errors (mainly, statistical) of these KMC results. 
Fig. fTTI shows that for the SF-1 state and relatively small 
precipitates, p < 10, such errors at the nucleation stage 
are significant, while at larger p, and at later stages of 
evolution, these errors decrease. For the LBS-1 model 
(Fig. [TO]) , differences between two KMC runs are lower 
as the fluctuation effects here are much weaker. Let us 
also note that the difference between values of d p at two 
neighboring curves, d Pl (t) and d P2 (t) in Figs. [TUTO! 
d(i) = (d pi — d P2 ), has the meaning of the density of 
precipitates that include i copper atoms with i between 
pi and p2- Pi < i < P2- 

Let us first discuss the KMCA results presented in 
Figs. [TT][TO1 First, we note that they qualitatively agree 
with the classical theory nucleation (CTN) ideas de- 
scribed in Sec. Ill A. In particular, all dependences d p (t) 
reveal presence of four main stages of decomposition: in- 
cubation, nucleation, growth and coarsening. For more 
quantitative comparison to the CTN, we should estimate 
the "critical" embryo size p c . In thermodynamic calcu- 
lations [27l - l29j . the critical embryo is precisely defined 
as the set of mean occupations {c{\ that corresponds to 
the saddle-point of the generalized free energy F{ci} in 
the multi-dimensional space a [28]. At the same time, 
for nucleating precipitates the analogous critical size can 
not be defined uniquely. It seems natural to define it 
as the lowest value of the embryo size p such that the 
most probable evolution path at p > p c is growth rather 
than shrinkage or dissolution, but due to the presence 
of fluctuations inherent to the nucleation process, such 
"kinetic" critical sizes can be determined just approxi- 
mately. 

Using KMCA results presented in Figs. [TTifMl we can 
estimate these p c from the shape of curves d p (t) at dif- 
ferent p. For the SF model, for which fluctuations in de- 
pendences d p (t) are pronounced at small p and decrease 
at larger p, we can suggest that beginning of decreasing 
these fluctuations with increasing p should correspond 
to the relation p > p c . For the states SF-1, SF-2 and 
SF-3, it corresponds to p c ~ 15, p c ~ 11, and p c ~ 8. 



FIG. 14: (color online). Same as in Fig. [TT]but for the LBS- 
1 state and p — 6, 11, 15 and 26, respectively, while dashed 
curve shows the SSA results for p = 15, go = 1.9. 
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As these estimates agree well with the "thermodynamic" 
critical sizes N c in Table II, we will use for the critical 
size p c in the SF-1, SF-2 and SF-3 alloy state the value 
15, 11 and 8, respectively. Note that in the SSA simula- 
tions, the p c values are used only to define the "end of 
nucleation" time t^f in Eq. ([7Tj) as the time of creation 
of the last critical precipitate, and slight variations of p c 
have almost no effect on evolution. 

Figs. ITlTfSI show that temporal dependences d p {t) at 
these p=p c have, generally, a common form similar to 
that obtained by SM [4J for their simplified model that 
agrees with the CTN ideas described in Sec. Ill A. At 
the same time, these dependences reveal at least two dif- 
ferences from the CTN ideas. 

(A) Fluctuations in dependences d Pc (t) are rather pro- 
nounced and significant even for "supercritical" embryos 
with p > p c , and the more so for undercritical ones with 
p < p c . It is particularly clear seen for the KMC results 
presented in Fig. [TT] by solid lines, while in Figs. [T2l 
and [TBI these fluctuations are partly smoothed due to the 
plotting of lesser number of the results. 

(B) These fluctuations are large and important not 
only during nucleation but also after its completion when 
the average density of supercritical precipitates ceases to 
increase. It is seen, in particular, in Fig. 1131 where curve 
d p (t) at p = 16 ~ 2p c reveals a wide and pronounced 
minimum between t = 3 • 10 5 and 4 • 10 5 sec falling here 
by about 1.5 times. Therefore, for the SF model, the fluc- 
tuations are important for evolution even at the growth 
stage. 

Let us now discuss dependences d p (t) for the LBS-1 
state for which the thermodynamic nucleation barrier F c 
(given in Table II) is low: F c < T. In this case, evolution 
of not large "thermodynamically supercritical" embryos 
with p > N c again differs from the CTN ideas: effects 
of fluctuations here are strong due to the low thermo- 
dynamic gain under growth of such embryos. As their 
stability is low, they often dissolve rather than steadily 
grow as the CTN supposes. It can qualitatively explain 
a significant decrease of "plateau" in curves d p {t) with 
p > 15 as compared to p=6 and p=10 in Fig. [TU For 
this case, it seems more adequate to suppose the kinetic 
critical size p c to significantly exceed N c , so that the 
dominant evolution path at p > p c would be growth of 
the embryo. For the LBS-1 state, we estimate this size 
as p c ~ 15. 

The SSA results for the precipitate density d p (t) pre- 
sented in Figs. ITTirnn have been obtained as described in 
Sec. IV. Note that rescaling of time t r used in the SSA 
to the time t used in the KMCA (illustrated by Fig. IT0|) 
corresponds to the above-described two-parametric fit of 
only "horizontal" intervals in Figs. [TTIfT4l while "ver- 
tical" intervals, that is, the density d p values, are cal- 
culated in the SSA with no adjustable parameters. We 
see that for all four alloy states considered, these d p (t) 
agree with the KMCA results within errors of these re- 
sults mentioned above. Fig. fTS] (as well as frame 115b 
below) also shows that changes in the SSA results due 



FIG. 15: (color online). Average number of copper atoms 
within a precipitate, (i), versus the evolution time t (in sec- 
onds). Frames a, b, c and d correspond to the alloy states 
SF-1, SF-2, SF-3 and LBS-1, respectively. Solid lines and 
squares correspond to the KMCA. Dashed lines correspond 
to the SSA with go from (|72[) . and chained line (in frame b), 
to go— 1.8. Arrows show the values of t that correspond to 
t? in Eq. CD}. 



FIG. 16: (color online). Upper frame: Distribution of copper 
atoms for the SF-1 alloy state just before nucleation, t~55 
sec, observed in the KMCA. Each sphere (blue online) corre- 
sponds to a copper atom. Lower frame: Distribution of local 
concentrations a for the same state as in the upper frame 
observed in the SSA. Relation between coloring and a value 
on the lattice site i is shown in the right part of the frame. 
Sites with a < 0.05 are not shown. 

to possible variations of go between go=1.7 and go=1.8 
mentioned in derivations of estimates (|72[) are relatively 
small. 

In Fig. [T5] we show temporal dependences of average 
number of copper atoms within a precipitate ("precip- 
itate size") (i)=i p (t). Results of the KMCA and SSA 
calculations seem usually to agree within the KMCA er- 
rors except for time intervals just after where i p SA 
notably exceed i p <MCA . This disagreement seems to be 
again related to crudeness of the oversimplified model 
(|7T|) for the length g(t r ) determining scale of fluctua- 
tive terms <5n/ in the stochastic equation (|56p . As men- 
tioned, model (I7T1) corresponds to the constant g(t r ) = 
go at t r < and to a sharp increase of g (and thus to 
practically abrupt switching-off fluctuations) at t r > t^f , 
while actually the length g starts to increase (and fluc- 
tuations Srif to decrease) immediately after beginning 
of nucleation, and g remains to be finite (and nucle- 
ation effects to be noticeable) at the growth stage, too. 
Thus at t r > t™, the SSA-calculated i p (t) in Fig. [Pol 
grow more rapidly than in the KMCA (and in reality) 
as actually this growth is still hampered by noticeable 
fluctuation effects. Therefore, we can expect that these 
disagreements will be reduced or vanished if the models 
of g(t r ) more realistic than (fTTj) will be used in the SSA. 

B. Features of microstructure at different stages of 
precipitation 

In Figs. [TMP?1 we show microstructure of the SF-1 al- 
loy state at different stages of evolution observed in the 
KMCA and SSA-bascd simulations. As the SSA simula- 



FIG. 17: (color online). Same as in Fig. [TS]but at the end of 
nucleation, t=1000 sec. Light spheres (yellow online) in the 
upper frame show copper atoms that belong to some super- 
critical precipitate. 
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FIG. 18: (color online). Same as in Fig. [17] but at the begin- 
ning of coarsening, £=2000 sec. 



FIG. 19: (color online). Same as in Fig. [17] but at some 
intermediate stage of coarsening, i=4000 sec. 

tion box was twice as large as that used in the KMCA, 
the total number of precipitates in the SS A-based (lower) 
frame of each of Figs. [TMT^l exceeds that presented in the 
KMC-based (upper) frame by about twice. Figs. [TB1[T51 
illustrate processes of nucleation, growth and coarsening 
discussed above. 

In comparison of the KMCA and SSA results in Figs. 
fTBlTT^l we note that the SSA, as well as any other statis- 
tical description, disregards details of particular atomic 
configurations presented in the KMCA snapshots, but it 
enables us to pick out essential features of microstructure 
which often can not be easily apprehended in these snap- 
shots. It is illustrated by Fig. [TBI where the SSA frame 
clearly shows a number of regions significantly enriched 
by copper atoms, that is, precursors of precipitates to 
be created, while such regions are not clearly seen in the 
analogous KMC frame. Similar differences in describing 
precipitates can be seen in Figs. [T7l[T9l Main quantita- 
tive characteristics of precipitates, such as their density 
and sizes discussed above, can be easily determined in 
both approaches. At the same time, the presence of sig- 
nificant crystalline anisotropy, particularly for not large 
precipitates, in the KMCA can be established only after 
some special analysis [TH, [l7|, while in the SSA frames 
it is seen at first sight. Difference between two descrip- 
tions is also evident for interfaces of precipitates. In the 
KMCA frames of Figs. [T?1[TT)I these interfaces are usually 
rather sharp but typically not flat and not regular, while 
in the analogous SSA frames, these interfaces seem to 
be somewhat diffuse and have typically "intermediate" 
values of local concentration: Cj ~ 0.5 (green online). 
The differences arise because each local concentration Cj 
in the SSA is obtained by averaging over some locally 
equilibrated vicinity of site i, that is, over a relatively 
rapid motion of surface atoms on the "not-filled" facet 
considered. At the same time, inner parts of precipi- 
tates (for both the KMCA and SSA) include only copper 
atoms which in the SSA frames is clearly seen on cuts of 
precipitates by the boundary planes of simulation box. 

Therefore, Figs. [T6HT91 also illustrate a complementary 
character of describing evolution of microstructure by the 
KMCA and the SSA. 



C. Kinetics of nucleation 

As discussed above, the SSA provides a partly averaged 
description of atomic distributions aiming mainly at ad- 
equate calculations of locally averaged quantities, such 
as the density of different precipitates, their structure, 



FIG. 20: (color online). Distribution of local concentrations 
within a plane containing several nucleating precipitates for 
the SF-1 state at the following reduced time t r values: (a) 
7.5, (b) 7.9, (c) 8.0, (d) 6.1, (e) 8.5, and (f) 9.0. For each 
point r of the figure, the c(r) value is obtained by inter- 
polation between a on neighboring lattice sites. Relation 
between coloring and a values is shown in the right part of 
the figure. 



morphology, etc. Therefore, quantitative treatments of 
phenomena which are mainly determined by fluctuations, 
such as processes of creation and evolution of under- 
critical and near-critical embryos, lies, generally, outside 
the scope of the SSA. However, as mentioned above, the 
qualitative changes of microstructure, such as those cor- 
responding to the nucleation process, can often be more 
easily followed in the SSA rather than in the KMCA. 
Therefore, it can be instructive to study kinetic details 
of nucleation with the use of the SSA, even though the 
scale of fluctuation effects in such study is most probably 
underestimated . 

Some results of such studies are presented in Figs. 
|2"0][2"21 which illustrate processes of sequent creation of 
three supercritical embryos for the SF-1 alloy state. 
Frames [20k-[20e also show processes of creation and dis- 
solution of an "undercritical" embryo: the concentration 
fluctuation in the right central part of these frames first 
increases up to the values Cj ~ 0.1, but then decreases 
and disappears. On the contrary, frames [20H-[20f illus- 
trate a "successful" nucleation process. The local fluctu- 
ation below that discussed above first increases in both 
size and amplitude, and then suddenly shrinks with for- 
mation of a "kinetic" supercritical embryo. Later on this 
embryo survives and grows, but this growth is first non- 
monotonic and includes a "partial dissolution" process 
illustrated by frames l2"Tk and l2"Tb . Note that at first 
stages of this process shown in frames [20b and [20f . the 
embryo is extended and shapeless, while later on it is 
rather wrong-shaped, thus it seems to have little in com- 
mon with the "thermodynamic" critical embryos shown 
in Fig. 

Creation and evolution of two other embryos shown in 
Figs. l2"Tlandl2"2"lproceeds similarly. In both cases, the ini- 
tial extended and shapeless fluctuation of concentration 
first suddenly shrinks so that maximum concentrations 
within it reach "critical" values Cj > 0.12, after which 
embryos seem to become "supercritical" . Then subse- 
quent fluctuations lead to a partial dissolution of these 
embryos, but in both cases they survive and later on start 
to grow. This growth correspond to "sucking" of copper 
atoms from adjacent regions and thus to depletion of cop- 
per concentration in these regions. 

The examples considered can also illustrate the kinetic 
mechanism of nucleation. It can be viewed as a local 
spinodal decomposition that starts when the amplitudes 
and extension of local fluctuative enrichment of concen- 
tration become large enough in order that the "uphill 
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FIG. 21: (color online). Same as in Fig. [20] but at the fol- 
lowing t r : (a) 21, (b) 24, (c) 29, (d) 31, (e) 34, and (f) 38. 



FIG. 22: (color online). Same as in Fig. [20] but at the fol- 
lowing t r : (a) 46, (b) 47, (c) 48, (d) 49, (e) 50, and (f) 53. 



diffusion" mechanism characteristic of spinodal decom- 
position [32I ] becomes operative. For the examples con- 
sidered, it seems to correspond to local concentrations 
Ci > 0.07-0.09 (while the uniform spinodal decomposition 
boundary, according to FigQ] is c s ~ 0.045) in the re- 
gion I > (7 — 8) a. Such interpretation of nucleation as a 
fluctuation-induced local spinodal decomposition can be 
useful for qualitative understanding of many phenomena 
in this field, in particular, of a great excess in probabili- 
ties of nucleation near the binodal curve observed exper- 
imentally [36| with respect to estimates of the classical 
theory of nucleation (|55|) . 

D. Changes in microstructural evolution under 
variations of temperature or concentration 

As mentioned in Sec. IV A, the kinetic type of al- 
loy decomposition is mainly determined by the value of 
the reduced supersaturation s (|65|) , which for metastable 
states under consideration is less than unity. Low s cor- 
respond to the states with the concentration and temper- 
ature (c, T) values close to the binodal curve for which 
a "deep NG" type of evolution with a low density and 
large sizes of nucleating precipitates is characteristic. An 
increase of s corresponds to approaching the spinodal 
curve and decreasing nucleation barriers, thus nucleating 
precipitates should become smaller (which is illustrated 
by Table II), while their density should increase. Basing 
on these considerations, we can expect that the reduced 
supersaturation s is the main parameter determining mi- 
crostructural evolution, but at the given s, microstruc- 
ture can also significantly vary with the concentration or 
temperature. 

To get an idea about these varuations, in Fig. [23] we 
show microstructure at the end of nucleation for the SF- 
1, SF-3 and SF-4 alloy states described by Tables I, II 
and Fig. [5] In simulations for the SF-4 state (which has 
has the same supersaturation as the SF-3 state but higher 
concentration and temperature) we used go — 1.9 value, 
same as for the SF-3 state. In accordance with consider- 
ations mentioned above, the increase of supersaturation 
s in states SF-3 and SF-4 with respect to SF-1 by about 



FIG. 23: (color online). Distribution of concentrations d at 
the end of nucleation for the following alloy states: (a) SF-1, 
(b) SF-3, and (c) SF-4. 



1.5 times leads to much higher density of nucleated pre- 
cipitates: by about 3.6 times for the SF-3 state, and by 
2.2 times for the SF-4 state. At the same time, frames 
l2"5b and l2l!b illustrate differences in microstructure for 
the same s but different c and T. For the state SF-4, 
precipitates are notably larger while their density is by 
1.5 times lower than those for the state SF-3. These dif- 
ferences can be qualitatively explained by the differences 
in characteristics of thermodynamic critical embryos for 
these two states presented in Table II and Fig. [5] both 
critical sizes and reduced nucleation barrier F c /T and for 
the SF-4 state notably exceed those for the SF-3 state. 



VI. CONCLUSIONS 

To conclude, we summarize the main results of this 
work. 

1. The consistent and computationally efficient 
stochastic statistical approach is developed to microscop- 
ically study kinetics of decomposition of metastable al- 
loys. 

2. In this approach, description of evolution in terms of 
certain reduced time includes no adjustable parameters. 
Rescaling of this reduced time to the physical time can 
usually be made with the use of few constants which can 
be estimated either from comparison to kinetic Monte 
Carlo simulations or from experiments. 

3. For several realistic models of iron-copper alloys 
studied, the results of this approach usually agree with 
the kinetic Monte Carlo simulation results within errors 
of these simulations. 

4. Oversimplified model (fTTj) for the important kinetic 
parameter of the theory, size of locally equilibrated re- 
gions, seems to be sufficient for describing main char- 
acteristics of microstructure. However, for an adequate 
description of temporal dependences we should use more 
realistic models discussed in Sees. IV B and V A. 
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